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We carry out a finite size scaling analysis of the jamming transition in frictionless bi-disperse 
soft core disks in two dimensions. We consider two different jamming protocols: (i) quench from 
random initial positions, and (ii) quasistatic shearing. By considering the fraction of jammed states 
as a function of packing fraction for systems with different numbers of particles, we determine the 
spatial correlation length critical exponent v w 1, and show that corrections to scaling are crucial for 
analyzing the data. We show that earlier numerical results yielding v < 1 are due to the improper 
neglect of these corrections. 
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Glassy behavior in condensed matter and granular sys- 
tems remains a topic of considerable controversy. In this 
context, the jamming of hard or soft core particles at zero 
temperature has been the focus of much recent effort. As 
the packing fraction <fr of a granular material increases, 
the system undergoes a sharp jamming transition from a 
fluid- like state to a rigid but disordered solid state [T]. It 
has been proposed that this T = transition is described 
by a critical point, with scaling behavior similar to that 
at a continuous phase transition as found in equilibrium 
systems [2 . A key signature of a continuous transition is 
a correlation length £ that diverges at the jamming <j)j, 
£ ~ — (j)j\~ u . Determination of the critical exponent v 
is thus a key goal in establishing and characterizing the 
critical nature of the jamming transition. 

While it has been suggested that the value of v is inde- 
pendent of the dimensionality of the system, or the spe- 
cific force law between particles [3] , the precise numerical 
value of v varies widely throughout the literature. From 
theoretical consideration of soft vibrational modes in the 
jammed solid, Wyart et al. [I] argued for v = 1/2. Nu- 
merical simulations of vibrational modes led Silbert et 
al. in two (2D) and three (3D) dimensions [5] to pos- 
tulate diverging transverse and longitudinal correlation 
lengths with exponents vt ~ 0.24 and Vt, ~ 0.48 respec- 
tively, fc-core percolation models, in mean field theory, 
also yield |6] two exponents v* — 1/4 and v# = 1/2, 
while a field theoretic approach [7 to jamming in 2D 
gave v = 1/4. Simulations by Drocco et al. [8] of a trace 
particle dragged through an incipient 2D jammed liquid 
resulted in a value v = 0.71±0.12, while from a numerical 
finite size scaling analysis of mechanically stable states in 
2D and 3D O'Hern et al. [3] found v = 0.71 ± 0.08. A 
scaling analysis of velocity correlations in simulated 2D 
shear driven flow by two of us [9| previously reported that 
v = 0.6 ±0.1. Hatano [10] obtained v = 0.73 ±0.05 from 
simulations of shear relaxation in 3D. while relaxation of 



random initial states to mechanical equilibrium in 2D led 
Head [II] to v = 0.57 ± 0.05. Heussinger and Barat [H] 
estimate v = 0.8 — 1.0 from displacement correlations in 
a 2D system under quasistatic shearing, while Heussinger 
et al. [13 find a dynamic correlation length in 2D with 
exponent v = 0.9. Establishing the precise value of v 
and determining whether all these correlations lengths 
are the same thus remains a crucial theoretical objective. 

In this work we present a detailed finite-size-scaling 
analysis of the jamming transition in frictionless bi- 
disperse soft core disks in 2D. Only through such a scaling 
analysis can one hope to clearly establish the singular be- 
havior of the system in the limit of infinite size, and the 
value of critical exponents. An advantage of the finitc- 
size-scaling method is that it allows one to compute the 
exponent v of the most divergent length scale without 
the need to explicitly calculate the correlation length £ 
itself. 

We consider two different jamming ensembles: (i) 
quench from random initial positions (RAND), and (ii) 
quasistatic shearing (QS) [H]. By considering the frac- 
tion of jammed states / as a function of packing fraction 
4> for systems with different numbers of particles N, we 
demonstrate that the correlation length critical exponent 
in both ensembles is v w 1 . We further show that correc- 
tions to scaling are crucial for understanding our data, 
and argue that earlier numerical results yielding v ss 0.7 
are due to the improper neglect of these corrections. Our 
results suggest that corrections to scaling may be im- 
portant in other scaling analyses of critical behavior at 
jamming, for example in rheological behavior. 

Our model is a 50:50 bi-disperse mixture of disks with 
diameters in the ratio 1.4 |3J. Particles interact with a 
soft core harmonic repulsion, 
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where r,-j is the distance between the centers of two par- 
ticles i and j, and dij is the sum of their radii. Length 
is in units such that the smaller diameter is unity, and 
energy is in units such that e = 1. A system of N disks 
in an area A thus has a packing fraction (density) 

= Att(0.5 2 + 0.7 2 )/(2A) . (2) 

We define our two ensembles as follow, (i) RAND: 
This is the ensemble introduced by O'Hern et al. [3]. 
We start with a fixed number of particles, N, at density 
(j), in a square box with periodic boundary conditions. 
Particles are put at random initial positions, and then a 
conjugate gradient method is used to relax the system to 
the nearest local energy minimum. The minima resulting 
from many such initial configurations (we use typically 
10000 for each value of 0) defines the ensemble, (ii) QS: 
At a fixed N and <f>, we start the system in a random 
initial configuration, and then apply a small shear strain 
step A7 using Lees-Edwards boundary conditions [14] , 
A conjugate gradient method then relaxes the system to 
the nearest local energy minimum, before the system is 
strained again by A7. The set of states obtained after the 
energy minimization, after a long total strain 7, defines 
the ensemble. We choose the strain step small enough 
that our results do not depend on the value of A7. For 
our biggest systems we use A7 = 10~ 5 . We average over 
10 — 20 independent runs, each sheared a total strain 
7^4 — 8; for our smaller sizes, we use 7 up to 200. 

In both ensembles, we stop the energy minimization 
when one of the following conditions is met (i) the rela- 
tive decrease in the energy after 50 iterations is smaller 
than 10 -10 , or (ii) the average energy per particle is 
E/N < 10~ 16 . In the latter case, we consider the re- 
sulting configuration to be unjammed. The key quantity 
in our analysis will be the fraction of jammed states in 
the ensemble at a given value of density, f((f>). We have 
verified that the energy bound (ii) gives a clear separa- 
tion between the jammed and unjammed states up to 
the largest system size we have studied. Further details 
of our numerical procedures may be found in Ref. |15j . 

In Fig. [T] we present our results for f((f>) for systems 
of varying number of particles N for both RAND and 
QS. We see that /(</>) sharpens up and approaches a step 
function in the limit N — > 00; this singularity in /(</>) as 
N — > 00 is characteristic of a quantity that has scaling 
dimension zero. We would thus expect, to leading order, 
the finite-size-scaling behavior, 

f(4>, L) = T (^L 1/l/ ) where 6<j> = <j> - fa , (3) 

<j)j is the jamming density in the thermodynamic limit 
N — > 00, v is the correlation length critical exponent, 
and L = y/~N is a measure of the linear size of the sys- 
tem. A key prediction of Eq. ^ is that at = fa, curves 
of f(4>,L) for different L should all intersect, having the 
common value .F(O); plotting f(<f),L) vs SifiL 1 /", curves 




FIG. 1: (color online) Fraction of jammed states / vs packing 
fraction (j>, for systems with number of particles N. (a) is the 
RAND ensemble, (b) is the QS ensemble. Insets show a blow 
up of the region where curves for different JV intersect. 



of different L should collapse to a common scaling curve. 
However careful inspection of our results in Fig. [I] (see in- 
sets) show that there is no common intersection point for 
the f((f>,L). This observation leads us to conclude that, 
for the sizes studied here, corrections to scaling must be 
included in our analysis. 

We can include such corrections to scaling by general- 
izing Eq. j3} to, 

f(<l>,L)=T (5<t>L 1 ^)+L- u Ti(s t l>L 1 ^ . (4) 

In the renormalization group framework for equilibrium 
critical phenomena, such corrections to scaling arise from 
a Taylor series expansion of the free energy in the leading 
irrelevant scaling field, whose scaling dimension is — uj 
[16] . We will define f c = ^(O) as the critical value of 
the jamming fraction at <pj in the limit L — > 00. 

One of the consequences of Eq. Q is that the func- 
tions /(</>, L) approach the L — > 00 limiting step function 
at different rates, depending on the value of /. If we 
define fa(L) as the value of 4> where f(<j),L) — f, then 
sufficiently close to fa we can expand the scaling func- 



FIG. 2: (color online) 4>f(L) vs L for different values of / for 
(a) RAND and (b) QS. Values of / increase from bottom to 
top. 



tions in Eq. Q to linear order in S(j) to obtain, 

cj )f (L) = ( f ) j-L- 1 ^[c () 5f-(c 1 -c 2 6f)L-"] , (5) 

where cq, c\, C2 are constants and Sf = f — f c . 

In Fig. [2] we plot 4>f(L) vs L for several values of /. 
To interpolate between our data points so as to define 
the values 4>}{L), we use the following procedure. We 
transform to a new variable F = ln[//(l — /)] and fit 
F((j)) to a fifth order polynomial over the range \F\ < 5. 
The result gives the solid lines in Fig. [T] We see that 
as / increases, 4>f{L) becomes non-monotonic, a clear 
signature of the change in sign of the leading term L^ 1 ^ 
111 Eq. as / increases above f c . We see that 4>j ~ 
0.8415 for RAND, while <j>j w 0.843 for QS. 

We consider next a determination of the exponent v 
via Eq. ([5| . To eliminate the imprecisely known value of 
</>,/, and to reduce the contribution from the correction 
to scaling given by C\, we consider the difference, 

w(L) = <f )f2 (L)~cl )fl (L)=aL- 1 ^(l + bL-") , (6) 

where both a and b are proportional to fa—fi. We choose 
fx and /2 symmetrically about f c (with f c as determined 
below), and plot w(L) vs L for RAND and QS in Fig. [3] 
We expect the correction term ~ to get smaller, and 
become negligible, as L increases. We therefore ignore 
the correction term and fit the data tow~ L~ x l v to get 
the solid line in Fig. [3| The insets show the resulting 
value of l/i/ as we drop successively smaller system sizes 
from the fit, fitting systems of size Af m in to N max (A^ max = 
16384 for RAND, N max = 4096 for QS, L = VN). As 
expected, the value of l/v saturates to a constant as N m - in 
increases and the effects of the correction term become 
negligible. We find from these fits the values l/v = 0.93± 
0.02 for RAND, and l/v = 0.91 ±0.02 for QS. If we then 
fit the data for all sizes to the full Eq. ([6]), including 
the correction term, we get values of l/v consistent with 
those above, however the estimated error in lo is too large 
to determine u> to any accuracy. 

To determine lo, and get a more accurate value for 
4>j, we use the following procedure. We fit the re- 
sults for 4>f{L) of Fig. [2] to a single power law <f>f(L) = 



FIG. 3: (color online) Width w(L) = (f>f 2 (L)-cj> fl (L) vs L for 
(a) RAND and (b) QS. We choose /1 and symmetrically 
about f c ; for RAND /1 = 0.7, f 2 = 0.9; for QS /1 = 0.5, 
/2 = 0.7. Straight line is a fit to w ~ L _1//l/ including all 
data. Insets show the fitted value of l/v as the minimum size 
system included in the fit, N m i U , is varied. 



cf)j — cLT r ' Vetf . Since 4>f{L) has such a single power law 
behavior only at f c , we expect that the % 2 of the fit will 
be smallest when / — f c . The fitted parameters at this 
f c then determine <pj and the exponent l/y e s = l/v + to. 
We show the results of such fits in Fig. |4j where we 
fit to system sizes Af m i n to N max , for the four differ- 
ent cases N m i n = 48, 64, 96, 128. We see that as A^ min 
increases, the low-/ side of the minimum in x 2 S e ts in- 
creasingly shallow. This is not surprising since the size of 
the correction term, and hence its effect on the fits, gets 
progressively smaller as N increases. Nevertheless we 
find quite stable values of the fitted parameters as Af m ; n 
varies, we find f c = 0.78±0.02, <j>j = 0.84177 ± 0.00001, 
l/v + uj = 1.7 ± 0.1 for RAND, and f c = 0.60 ± 0.03, 
4>j = 0.8432±0.0001, 1/v + uj = 1.85±0.03 for QS. Com- 
bining with our earlier results for l/v we get uj = 0.8 ±0.1 
for RAND and w = 0.94 ± 0.05 for QS. The solid lines 
in Fig. [2] result from fits to Eq. ([5| where we have fixed 
1/u and l/v + lo to the values found from the analyses of 
Figs. [3] and gj 

It is interesting to compare our results against the fi- 
nite size scaling analysis of O'Hern et al. [3], who con- 
sidered the RAND ensemble. In that work, the authors 
considered the distribution P(<fi,L) — df(<fi,L)/d<p, the 
probability density for a system of size L to have its 
particular jamming density at <fi. By considering how 
the location 4>o(L) of the peak in P(<f>,L) approached its 
L — > 00 limit <pj, the authors defined the critical expo- 
nent 'V by, 4>j — 4>q ~ L _1 /' v ", and found the value 
'V = 0.71 ± 0.08. In terms of our analysis, we see that 
4>o is the same as our <pf , where /o locates the steep- 
est slope of /(</>, L), and V is just our u e s. In light of 
corrections to scaling, we see that V" should not be iden- 
tified as the correlation length exponent; it is an effective 
exponent that arises from fitting (j) to single power law, 
when the true behavior as in Eq. |5]) is governed by two 
different power laws with exponents l/v and l/v + lo. If 
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FIG. 4: (color online) Results from fitting data of Fig. [5] to 

4>f{L) — <j)j — cL _1/ ' 1 ' cff , using system sizes A m i n to iV max . 
We show results for jV min = 48, 64, 96, f 28. Panels (a)-(c) are 
for RAND, panels (d)-(f) are for QS. (a), (d) is the X 2 /dof of 
the fit (dof = number of data points minus number of fitting 
parameters); the minimum of x 2 /dof locates f c , which then 
determines cj>j, as shown in (b), (e), and the value of l/v e s = 
l/v + uj, as shown in (c), (f). Insets show the dependence of 
the fitted parameters on the value of -/V m i n . 



we take /o = 0.5, our Fig. [if: for the case iV m ; n = 64 gives 
l/^cff = 1-32, or z/ ff = 0.76, in good agreement with the 
value found by O'Hern et al. 

O'Hern et al. similarly define the full width at half 
maximum of P((j),L), w(L), and find the scaling w ~ 
N' n ~ L- 2Q , with fi = 0.55 ± 0.03 or 2ft = 1.10 ± 0.06. 
With suitable choices of /i and fi , this w is the same as 
our w of Eq. and hence we expect for asymptotically 
large N (where corrections to scaling become negligible) 
to find 2fi = \jv. If we assume a Gaussian form for 
P(4>, L) then the full width criterion corresponds to /i = 
0.124 and — 0.876. Computing this w and fitting using 
sizes N = 64 to 4096, the same range as O'Hern et al., we 
get 2f2 = 1.035 ± 0.002, in agreement with O'Hern et al. 
within their estimated errors. However if we use up to 
our largest size N ~ 16384, then increase iV m i n , we find 
f2 to systematically decrease, becoming 2f2 = 0.97 ± 0.01 
when iV m i n = 2048. Our result remains larger than the 
1/v = 0.93 found in Fig. [3j perhaps because f\ is so far 
from f c that additional corrections to scaling arise from 
non-linearities in the scaling functions. Thus we conclude 
that it is O'Hern et al.'s 1 / (20) that is asymptotically the 
correlation length exponent v, rather than their 'V (our 
24fif)) and that their value for 2f2 is larger than our l/i> 
due to their more limited range of sizes and their neglect 
of corrections to scaling. 

In other recent work [17], it is found that corrections 
to scaling must similarly be included to properly describe 
the critical scaling of rheology under applied shear strain 
rate 7, for rates of the size typically used in simulations. 
The value <j)j = 0.8415 for shear driven jamming that was 



reported in earlier work by two of us [9] , is lower than the 
corresponding <j)j = 0.8432 found for QS here, due the 
neglect in that work of corrections to scaling. Similarly, 
the low value v — 0.6 reported in that work also results 
from the earlier failure to include corrections to scaling. 
We expect that other numerically reported values of <j)j 
and v may similarly be inaccurate due to the neglect of 
corrections to scaling in the analysis. 

To conclude, we have demonstrated that, for the sizes 
N typically used in simulations, including corrections to 
scaling is crucial for a proper description of the critical 
behavior at jamming, in both RAND and QS ensembles. 
Although we know no apriori reason why this should be 
so, it is interesting to note that corrections to scaling 
are similarly important in spin glass problems, another 
system in which the "ordered" state appears spatially 
random [16] . Within our estimated accuracy we find the 
correlation length exponent v and the correction to scal- 
ing exponent uj to be roughly equal for the two ensem- 
bles. We find uj = 0.89 ± 0.12, and 1/v = 0.92 ± 0.02, or 
v = 1.09 ± 0.02. While the estimated statistical error in 
v is small [18] . our range of system sizes L is not suffi- 
ciently large for us rule out the possibility that systematic 
errors, for example from additional or higher order cor- 
rections to scaling, could slightly alter the value of these 
exponents to v = oj = 1. 
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